#include "common.h"
#include "Altitude.h"
#include "AtmosModel.h"
#include "AerosolConcentration.h"

/* Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the target 
   is not at sea level.


	Given the altitude of the target in kilometers as input, we transform the
	original atmospheric profile (Pressure, Temperature, Water Vapor, Ozone) 
	so that first level of the new profile is the one at the target altitude. 
	We also compute the new integrated content in water vapor and ozone, that
	are used as outputs or in computations when the user chooses to enter a
	specific amount of Ozone and Water Vapor.
*/
void Altitude::pressure(AtmosModel& atms, float& uw, float& uo3)
{
	/* log linear interpolation */
	if(xps >= 100) xps = 99.99f;
		
	int i;
	for(i = 0; atms.z[i] <= xps; i++);
		
	int isup = i;
	int iinf = i - 1;

	float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
	float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
	float ps = (float)exp((xps - xb) / xa);

	/* interpolating temperature wator vapor and ozone profile versus altitude */
	float xalt = xps;
	float xtemp = (atms.t[isup] - atms.t[iinf]) / (atms.z[isup] - atms.z[iinf]);
	xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
	float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
	xwo = xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
	float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
	xwh = xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];

	/* updating atmospheric profile
	  1rst level: target     , complete to 34
	  with interpolated layers */
	atms.z[0] = xalt;                                                          
	atms.p[0] = ps;
	atms.t[0] = xtemp;
	atms.wh[0] = xwh;
	atms.wo[0] = xwo;

	for (i = 1; i < 33 - iinf; ++i)
    {
		atms.z[i] = atms.z[i + iinf];
		atms.p[i] = atms.p[i + iinf];
		atms.t[i] = atms.t[i + iinf];
		atms.wh[i] = atms.wh[i + iinf];
		atms.wo[i] = atms.wo[i + iinf];
    }
	int l = 33 - iinf - 1;
	for (i = l; i < 34; ++i)
    {
      atms.z[i] = (atms.z[33] - atms.z[l]) * (i - l) / (33 - l) + atms.z[l];
      atms.p[i] = (atms.p[33] - atms.p[l]) * (i - l) / (33 - l) + atms.p[l];
      atms.t[i] = (atms.t[33] - atms.t[l]) * (i - l) / (33 - l) + atms.t[l];
      atms.wh[i] = (atms.wh[33] - atms.wh[l]) * (i - l) / (33 - l) + atms.wh[l];
      atms.wo[i] = (atms.wo[33] - atms.wo[l]) * (i - l) / (33 - l) + atms.wo[l];
    }

	/* compute modified h2o and o3 integrated content */
	uw = 0;
	uo3 = 0;
	const float g = 98.1f;
	const float air = 0.028964f/0.0224f;
	const float ro3 = 0.048f/0.0224f;

	float rmwh[34];
	float rmo3[34];
	int k;
	for (k = 0; k < 33; ++k)
    {
		float roair = air * 273.16f * atms.p[k] / (atms.t[k] * 1013.25f);
		rmwh[k] = atms.wh[k] / (roair * 1e3f);
		rmo3[k] = atms.wo[k] / (roair * 1e3f);
    }

	for (k = 1; k < 33; ++k)
    {
		float ds = (atms.p[k - 1] - atms.p[k]) / atms.p[0];
		uw += (rmwh[k] + rmwh[k - 1]) * ds / 2.f;
		uo3 += (rmo3[k] + rmo3[k - 1]) * ds / 2.f;
    }
	uw = uw * atms.p[0] * 100.f / g;
	uo3 = uo3 * atms.p[0] * 100.f / g;
	uo3 = uo3 * 1e3f / ro3;
}

/*
Function: Update the atmospheric profile (P(z),T(z),H2O(z),O3(z)) in case the observer is on
board an aircraft.

Description: Given the altitude or pressure at aircraft level as input, the first task is to
compute the altitude (in case the pressure has been entered) or the pressure (in case the altitude has
been entered) at plane level. Then, a new atmospheric profile is created (Pp,Tp,H2Op,O3p) for which
the last level is located at the plane altitude. This profile is used in the gaseous absorption
computation (ABSTRA.f) for the path from target to sensor (upward transmission). The ozone and
water vapor integrated content of the "plane" atmospheric profile are also an output of this
subroutine. The last output is the proportion of molecules below plane level which is useful in
scattering computations (OS.f,ISO.f).
*/
void Altitude::presplane(AtmosModel& atms)
{
	/* log linear interpolation */
	xpp += atms.z[0];
	if(xpp >= 100) xpp = 1000;

	int i;
	for(i = 0; atms.z[i] <= xpp; i++);

	int isup = i;
	int iinf = i-1;

	float xa = (float)((atms.z[isup] - atms.z[iinf]) / log(atms.p[isup] / atms.p[iinf]));
	float xb = (float)(atms.z[isup] - xa * log(atms.p[isup]));
	float ps = (float)(exp((xpp - xb) / xa));

	/* interpolating temperature wator vapor and ozone profile versus altitud */
	float xalt = xpp;
	float xtemp  = (atms.t[isup] - atms.t[iinf])/ (atms.z[isup] - atms.z[iinf]);
	xtemp = xtemp * (xalt - atms.z[iinf]) + atms.t[iinf];
	float xwo = (atms.wo[isup] - atms.wo[iinf]) / (atms.z[isup] - atms.z[iinf]);
	xwo =  xwo * (xalt - atms.z[iinf]) + atms.wo[iinf];
	float xwh = (atms.wh[isup] - atms.wh[iinf]) / (atms.z[isup] - atms.z[iinf]);
	xwh =  xwh * (xalt - atms.z[iinf]) + atms.wh[iinf];

	/* updating atmospheric profile
	  last level: plane     , complete to 34
	  with interpolated layers */
	for(i = 0; i <= iinf; i++)
	{
		plane_sim.zpl[i] = atms.z[i];
		plane_sim.ppl[i] = atms.p[i];
		plane_sim.tpl[i] = atms.t[i];
		plane_sim.whpl[i] = atms.wh[i];
		plane_sim.wopl[i] = atms.wo[i];
	}

	for(i = iinf+1; i < 34; i++)
	{
		plane_sim.zpl[i] = xalt;
		plane_sim.ppl[i] = ps;
		plane_sim.tpl[i] = xtemp;
		plane_sim.whpl[i] = xwh;
		plane_sim.wopl[i] = xwo;
	}

	/* compute modified h2o and o3 integrated content
	 compute conversion factor for rayleigh optical thickness computation
	 ftray=rp/rt */
	atms.uw = 0;
	atms.uo3 = 0;
	const float g = 98.1f;
	const float air = 0.028964f/0.0224f;
	const float ro3 = 0.048f/0.0224f;
	float rt = 0;
	float rp = 0;

	float rmo3[34];
	float rmwh[34];
	int k;
	for(k = 0; k < 33; k++)
	{
		float roair = (float)(air * 273.16 * plane_sim.ppl[k] / (1013.25 * plane_sim.tpl[k]));
		rmwh[k] = atms.wh[k] / (roair * 1000);
		rmo3[k] = atms.wo[k] / (roair * 1000);
		rt += (atms.p[k+1] / atms.t[k+1] + atms.p[k] / atms.p[k]) * (atms.z[k+1] - atms.z[k]);
		rp += (plane_sim.ppl[k+1] / plane_sim.tpl[k+1] + plane_sim.ppl[k] / plane_sim.tpl[k]) 
			* (plane_sim.zpl[k+1] - plane_sim.zpl[k]);
	}

	ftray = rp / rt;
	for(k = 1; k < 33; k++)
	{
		float ds = (plane_sim.ppl[k-1] - plane_sim.ppl[k]) / plane_sim.ppl[0];
		atms.uw += (rmwh[k] + rmwh[k-1])*ds/2;
		atms.uo3+= (rmo3[k] + rmo3[k-1])*ds/2;
	}

	atms.uw *= plane_sim.ppl[0] * 100 / g;
	atms.uo3*= plane_sim.ppl[0] * 100 / g;
	atms.uo3*= 1000 / ro3;
}

void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
{
    xps = original_xps;
    xpp = original_xpp;

	float uwus;
	float uo3us;
	if(xps <= 0)
	{
		xps = 0;
		uwus = 1.424f;
		uo3us = 0.344f;
	}
	else if(atms.idatm != 8) pressure(atms, atms.uw, atms.uo3);
	else pressure(atms, uwus, uo3us);

	if(xpp <= 0)
	{
		/* ground measurement option */
		palt = 0;
		pps = atms.p[0];
		idatmp = 0;
		taer55p = 0;
		puw = 0;
	}
	else if(xpp >= 100)
	{
		/* satellite case of equivalent */
		palt = 1000;
		pps = 0;
		taer55p = aerocon.taer55;
		puw = 0;
		ftray = 1;
		idatmp = 4;
	}
	else
	{
		/* "real" plane case */
		cin >> puw;
		cin >> puo3;
		cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
		if ( puw < 0 )
		{
			presplane(atms);
			idatmp = 2;

			if (atms.idatm == 8)
			{
				puwus = puw;
				puo3us = puo3;
				puw *= atms.uw / uwus;
				puo3 *= atms.uo3 / uo3us;
				idatmp = 8;
			}
		}
		else
		{
			presplane(atms);
			idatmp = 8;
		}

		palt = plane_sim.zpl[33] - atms.z[0];
		pps = plane_sim.ppl[33];
		cin >> taer55p;

		if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03))
		{
			/* a scale heigh of 2km is assumed in case no value is given for taer55p */
			taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
		}
		else
		{
			/* compute effective scale heigh */
			double sham = exp(-palt / 4);
			double sha = 1 - (taer55p / aerocon.taer55);

			if( sha >= sham) taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
			else {
				sha = -palt/log(sha);
				taer55p = (float)(aerocon.taer55 * (1 - exp(-palt/sha)));
			}
		}
	}
}

void Altitude::parse()
{
	cin >> original_xps;
	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
    original_xps = -original_xps;
    
	cin >> original_xpp;
	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
    original_xpp = -original_xpp;
}

/* --- plane simulation output if selected ---- */
void Altitude::print()
{
	if(palt < 1000)
	{
		Output::Ln();
		Output::WriteLn(22," plane simulation description ");
		Output::WriteLn(22," ---------------------------- ");
		
		ostringstream s1;
		s1.setf(ios::fixed, ios::floatfield);
		s1.precision(2);
		s1 << " plane  pressure          [mb] " << setw(9) << pps << ends;
		Output::WriteLn(10,s1.str());

		ostringstream s2;
		s2.setf(ios::fixed, ios::floatfield);
		s2.precision(3);
		s2 << " plane  altitude absolute [km] " << setw(9) << plane_sim.zpl[33] << ends;
		Output::WriteLn(10,s2.str());

		
		Output::WriteLn(15," atmosphere under plane description: ");

		ostringstream s3;
		s3.setf(ios::fixed, ios::floatfield);
		s3.precision(3);
		s3 << " ozone content            " << setw(9) << puo3 << ends;
		Output::WriteLn(15,s3.str());


		ostringstream s4;
		s4.setf(ios::fixed, ios::floatfield);
		s4.precision(3);
		s4 << " h2o   content            " << setw(9) << puw << ends;
		Output::WriteLn(15,s4.str());

		ostringstream s5;
		s5.setf(ios::fixed, ios::floatfield);
		s5.precision(3);
		s5 << "aerosol opt. thick. 550nm " << setw(9) << taer55p << ends;
		Output::WriteLn(15,s5.str());
	}
}

Altitude Altitude::Parse()
{
	Altitude alt;
	alt.parse();
	return alt;
}
